//
// File: eig.cpp
//
// MATLAB Coder version            : 5.1
// C/C++ source code generated on  : 22-Feb-2021 03:11:30
//

// Include Files
#include "eig.h"
#include "FLAE_data.h"
#include "FLAE_rtwutil.h"
#include "rt_nonfinite.h"
#include "xhseqr.h"
#include "xnrm2.h"
#include "xzggev.h"
#include "xzlarf.h"
#include "rt_nonfinite.h"
#include <cmath>
#include <cstring>

// Function Definitions
//
// Arguments    : const double A[16]
//                creal_T V[16]
//                creal_T D[16]
// Return Type  : void
//
namespace coder
{
  void eig(const double A[16], creal_T V[16], creal_T D[16])
  {
    creal_T At[16];
    creal_T alpha1[4];
    creal_T beta1[4];
    double b_A[16];
    double b_V[16];
    double work[4];
    double tau[3];
    double colnorm;
    double scale;
    double t;
    int ix0;
    int k;
    int knt;
    boolean_T p;
    p = true;
    for (k = 0; k < 16; k++) {
      if ((!p) || (rtIsInf(A[k]) || rtIsNaN(A[k]))) {
        p = false;
      }
    }

    if (!p) {
      for (int i = 0; i < 16; i++) {
        V[i].re = rtNaN;
        V[i].im = 0.0;
        D[i].re = 0.0;
        D[i].im = 0.0;
      }

      D[0].re = rtNaN;
      D[0].im = 0.0;
      D[5].re = rtNaN;
      D[5].im = 0.0;
      D[10].re = rtNaN;
      D[10].im = 0.0;
      D[15].re = rtNaN;
      D[15].im = 0.0;
    } else {
      int b_i;
      int exitg1;
      int j;
      boolean_T exitg2;
      p = true;
      j = 0;
      exitg2 = false;
      while ((!exitg2) && (j < 4)) {
        b_i = 0;
        do {
          exitg1 = 0;
          if (b_i <= j) {
            if (!(A[b_i + (j << 2)] == A[j + (b_i << 2)])) {
              p = false;
              exitg1 = 1;
            } else {
              b_i++;
            }
          } else {
            j++;
            exitg1 = 2;
          }
        } while (exitg1 == 0);

        if (exitg1 == 1) {
          exitg2 = true;
        }
      }

      if (p) {
        int i;
        p = true;
        for (k = 0; k < 16; k++) {
          if ((!p) || (rtIsInf(A[k]) || rtIsNaN(A[k]))) {
            p = false;
          }
        }

        if (!p) {
          for (i = 0; i < 16; i++) {
            b_V[i] = rtNaN;
          }

          knt = 2;
          for (j = 0; j < 3; j++) {
            if (knt <= 4) {
              std::memset(&b_V[(j * 4 + knt) + -1], 0, (5 - knt) * sizeof(double));
            }

            knt++;
          }

          for (i = 0; i < 16; i++) {
            b_A[i] = rtNaN;
          }
        } else {
          int rowleft;
          std::memcpy(&b_A[0], &A[0], 16U * sizeof(double));
          work[0] = 0.0;
          work[1] = 0.0;
          work[2] = 0.0;
          work[3] = 0.0;
          for (b_i = 0; b_i < 3; b_i++) {
            double absxk;
            int alpha1_tmp_tmp;
            int im1n_tmp_tmp;
            int in;
            int lastc;
            int lastv;
            im1n_tmp_tmp = b_i << 2;
            in = (b_i + 1) << 2;
            alpha1_tmp_tmp = (b_i + im1n_tmp_tmp) + 1;
            t = b_A[alpha1_tmp_tmp];
            if (b_i + 3 < 4) {
              lastv = b_i + 1;
            } else {
              lastv = 2;
            }

            ix0 = (lastv + im1n_tmp_tmp) + 2;
            tau[b_i] = 0.0;
            scale = internal::blas::xnrm2(2 - b_i, b_A, ix0);
            if (scale != 0.0) {
              scale = rt_hypotd_snf(t, scale);
              if (t >= 0.0) {
                scale = -scale;
              }

              if (std::abs(scale) < 1.0020841800044864E-292) {
                knt = -1;
                i = (ix0 - b_i) + 1;
                do {
                  knt++;
                  for (k = ix0; k <= i; k++) {
                    b_A[k - 1] *= 9.9792015476736E+291;
                  }

                  scale *= 9.9792015476736E+291;
                  t *= 9.9792015476736E+291;
                } while (!(std::abs(scale) >= 1.0020841800044864E-292));

                scale = rt_hypotd_snf(t, internal::blas::xnrm2(2 - b_i, b_A, ix0));
                if (t >= 0.0) {
                  scale = -scale;
                }

                tau[b_i] = (scale - t) / scale;
                absxk = 1.0 / (t - scale);
                i = (ix0 - b_i) + 1;
                for (k = ix0; k <= i; k++) {
                  b_A[k - 1] *= absxk;
                }

                for (k = 0; k <= knt; k++) {
                  scale *= 1.0020841800044864E-292;
                }

                t = scale;
              } else {
                tau[b_i] = (scale - t) / scale;
                absxk = 1.0 / (t - scale);
                i = (ix0 - b_i) + 1;
                for (k = ix0; k <= i; k++) {
                  b_A[k - 1] *= absxk;
                }

                t = scale;
              }
            }

            b_A[alpha1_tmp_tmp] = 1.0;
            k = in + 1;
            if (tau[b_i] != 0.0) {
              lastv = 2 - b_i;
              knt = (alpha1_tmp_tmp - b_i) + 2;
              while ((lastv + 1 > 0) && (b_A[knt] == 0.0)) {
                lastv--;
                knt--;
              }

              lastc = 4;
              exitg2 = false;
              while ((!exitg2) && (lastc > 0)) {
                rowleft = in + lastc;
                ix0 = rowleft;
                do {
                  exitg1 = 0;
                  if (ix0 <= rowleft + (lastv << 2)) {
                    if (b_A[ix0 - 1] != 0.0) {
                      exitg1 = 1;
                    } else {
                      ix0 += 4;
                    }
                  } else {
                    lastc--;
                    exitg1 = 2;
                  }
                } while (exitg1 == 0);

                if (exitg1 == 1) {
                  exitg2 = true;
                }
              }
            } else {
              lastv = -1;
              lastc = 0;
            }

            if (lastv + 1 > 0) {
              int i1;
              int ix;
              if (lastc != 0) {
                if (0 <= lastc - 1) {
                  std::memset(&work[0], 0, lastc * sizeof(double));
                }

                ix = alpha1_tmp_tmp;
                i = (in + (lastv << 2)) + 1;
                for (knt = k; knt <= i; knt += 4) {
                  rowleft = 0;
                  i1 = (knt + lastc) - 1;
                  for (ix0 = knt; ix0 <= i1; ix0++) {
                    work[rowleft] += b_A[ix0 - 1] * b_A[ix];
                    rowleft++;
                  }

                  ix++;
                }
              }

              if (!(-tau[b_i] == 0.0)) {
                knt = in;
                rowleft = alpha1_tmp_tmp;
                for (j = 0; j <= lastv; j++) {
                  if (b_A[rowleft] != 0.0) {
                    absxk = b_A[rowleft] * -tau[b_i];
                    ix = 0;
                    i = knt + 1;
                    i1 = lastc + knt;
                    for (k = i; k <= i1; k++) {
                      b_A[k - 1] += work[ix] * absxk;
                      ix++;
                    }
                  }

                  rowleft++;
                  knt += 4;
                }
              }
            }

            internal::reflapack::xzlarf(3 - b_i, 3 - b_i, (b_i + im1n_tmp_tmp) +
              2, tau[b_i], b_A, (b_i + in) + 2, work);
            b_A[alpha1_tmp_tmp] = t;
          }

          std::memcpy(&b_V[0], &b_A[0], 16U * sizeof(double));
          for (j = 2; j >= 0; j--) {
            ix0 = (j + 1) << 2;
            for (b_i = 0; b_i <= j; b_i++) {
              b_V[ix0 + b_i] = 0.0;
            }

            i = j + 3;
            for (b_i = i; b_i < 5; b_i++) {
              knt = ix0 + b_i;
              b_V[knt - 1] = b_V[knt - 5];
            }
          }

          b_V[1] = 0.0;
          b_V[2] = 0.0;
          b_V[3] = 0.0;
          b_V[0] = 1.0;
          knt = 2;
          work[0] = 0.0;
          work[1] = 0.0;
          work[2] = 0.0;
          work[3] = 0.0;
          for (b_i = 2; b_i >= 0; b_i--) {
            rowleft = (b_i + (b_i << 2)) + 5;
            if (b_i + 1 < 3) {
              b_V[rowleft] = 1.0;
              internal::reflapack::xzlarf(3 - b_i, 2 - b_i, rowleft + 1, tau[knt],
                b_V, rowleft + 5, work);
              ix0 = rowleft + 2;
              i = (rowleft - b_i) + 3;
              for (k = ix0; k <= i; k++) {
                b_V[k - 1] *= -tau[knt];
              }
            }

            b_V[rowleft] = 1.0 - tau[knt];
            for (j = 0; j < b_i; j++) {
              b_V[(rowleft - j) - 1] = 0.0;
            }

            knt--;
          }

          internal::lapack::xhseqr(b_A, b_V);
        }

        for (j = 0; j < 3; j++) {
          b_A[(j + (j << 2)) + 1] = 0.0;
          for (b_i = 0; b_i <= j; b_i++) {
            b_A[b_i + ((j + 1) << 2)] = 0.0;
          }
        }

        for (i = 0; i < 16; i++) {
          V[i].re = b_V[i];
          V[i].im = 0.0;
          D[i].re = b_A[i];
          D[i].im = 0.0;
        }
      } else {
        double absxk;
        for (int i = 0; i < 16; i++) {
          At[i].re = A[i];
          At[i].im = 0.0;
        }

        internal::reflapack::xzggev(At, &knt, alpha1, beta1, V);
        colnorm = 0.0;
        scale = 3.3121686421112381E-170;
        for (k = 1; k < 5; k++) {
          absxk = std::abs(V[k - 1].re);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }

          absxk = std::abs(V[k - 1].im);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }
        }

        colnorm = scale * std::sqrt(colnorm);
        for (j = 1; j < 5; j++) {
          scale = V[j - 1].re;
          t = V[j - 1].im;
          if (t == 0.0) {
            absxk = scale / colnorm;
            scale = 0.0;
          } else if (scale == 0.0) {
            absxk = 0.0;
            scale = t / colnorm;
          } else {
            absxk = scale / colnorm;
            scale = t / colnorm;
          }

          V[j - 1].re = absxk;
          V[j - 1].im = scale;
        }

        colnorm = 0.0;
        scale = 3.3121686421112381E-170;
        for (k = 5; k < 9; k++) {
          absxk = std::abs(V[k - 1].re);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }

          absxk = std::abs(V[k - 1].im);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }
        }

        colnorm = scale * std::sqrt(colnorm);
        for (j = 5; j < 9; j++) {
          scale = V[j - 1].re;
          t = V[j - 1].im;
          if (t == 0.0) {
            absxk = scale / colnorm;
            scale = 0.0;
          } else if (scale == 0.0) {
            absxk = 0.0;
            scale = t / colnorm;
          } else {
            absxk = scale / colnorm;
            scale = t / colnorm;
          }

          V[j - 1].re = absxk;
          V[j - 1].im = scale;
        }

        colnorm = 0.0;
        scale = 3.3121686421112381E-170;
        for (k = 9; k < 13; k++) {
          absxk = std::abs(V[k - 1].re);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }

          absxk = std::abs(V[k - 1].im);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }
        }

        colnorm = scale * std::sqrt(colnorm);
        for (j = 9; j < 13; j++) {
          scale = V[j - 1].re;
          t = V[j - 1].im;
          if (t == 0.0) {
            absxk = scale / colnorm;
            scale = 0.0;
          } else if (scale == 0.0) {
            absxk = 0.0;
            scale = t / colnorm;
          } else {
            absxk = scale / colnorm;
            scale = t / colnorm;
          }

          V[j - 1].re = absxk;
          V[j - 1].im = scale;
        }

        colnorm = 0.0;
        scale = 3.3121686421112381E-170;
        for (k = 13; k < 17; k++) {
          absxk = std::abs(V[k - 1].re);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }

          absxk = std::abs(V[k - 1].im);
          if (absxk > scale) {
            t = scale / absxk;
            colnorm = colnorm * t * t + 1.0;
            scale = absxk;
          } else {
            t = absxk / scale;
            colnorm += t * t;
          }
        }

        colnorm = scale * std::sqrt(colnorm);
        for (j = 13; j < 17; j++) {
          scale = V[j - 1].re;
          t = V[j - 1].im;
          if (t == 0.0) {
            absxk = scale / colnorm;
            scale = 0.0;
          } else if (scale == 0.0) {
            absxk = 0.0;
            scale = t / colnorm;
          } else {
            absxk = scale / colnorm;
            scale = t / colnorm;
          }

          V[j - 1].re = absxk;
          V[j - 1].im = scale;
        }

        std::memset(&D[0], 0, 16U * sizeof(creal_T));
        if (beta1[0].im == 0.0) {
          if (alpha1[0].im == 0.0) {
            D[0].re = alpha1[0].re / beta1[0].re;
            D[0].im = 0.0;
          } else if (alpha1[0].re == 0.0) {
            D[0].re = 0.0;
            D[0].im = alpha1[0].im / beta1[0].re;
          } else {
            D[0].re = alpha1[0].re / beta1[0].re;
            D[0].im = alpha1[0].im / beta1[0].re;
          }
        } else if (beta1[0].re == 0.0) {
          if (alpha1[0].re == 0.0) {
            D[0].re = alpha1[0].im / beta1[0].im;
            D[0].im = 0.0;
          } else if (alpha1[0].im == 0.0) {
            D[0].re = 0.0;
            D[0].im = -(alpha1[0].re / beta1[0].im);
          } else {
            D[0].re = alpha1[0].im / beta1[0].im;
            D[0].im = -(alpha1[0].re / beta1[0].im);
          }
        } else {
          t = std::abs(beta1[0].re);
          scale = std::abs(beta1[0].im);
          if (t > scale) {
            scale = beta1[0].im / beta1[0].re;
            absxk = beta1[0].re + scale * beta1[0].im;
            D[0].re = (alpha1[0].re + scale * alpha1[0].im) / absxk;
            D[0].im = (alpha1[0].im - scale * alpha1[0].re) / absxk;
          } else if (scale == t) {
            if (beta1[0].re > 0.0) {
              scale = 0.5;
            } else {
              scale = -0.5;
            }

            if (beta1[0].im > 0.0) {
              absxk = 0.5;
            } else {
              absxk = -0.5;
            }

            D[0].re = (alpha1[0].re * scale + alpha1[0].im * absxk) / t;
            D[0].im = (alpha1[0].im * scale - alpha1[0].re * absxk) / t;
          } else {
            scale = beta1[0].re / beta1[0].im;
            absxk = beta1[0].im + scale * beta1[0].re;
            D[0].re = (scale * alpha1[0].re + alpha1[0].im) / absxk;
            D[0].im = (scale * alpha1[0].im - alpha1[0].re) / absxk;
          }
        }

        if (beta1[1].im == 0.0) {
          if (alpha1[1].im == 0.0) {
            D[5].re = alpha1[1].re / beta1[1].re;
            D[5].im = 0.0;
          } else if (alpha1[1].re == 0.0) {
            D[5].re = 0.0;
            D[5].im = alpha1[1].im / beta1[1].re;
          } else {
            D[5].re = alpha1[1].re / beta1[1].re;
            D[5].im = alpha1[1].im / beta1[1].re;
          }
        } else if (beta1[1].re == 0.0) {
          if (alpha1[1].re == 0.0) {
            D[5].re = alpha1[1].im / beta1[1].im;
            D[5].im = 0.0;
          } else if (alpha1[1].im == 0.0) {
            D[5].re = 0.0;
            D[5].im = -(alpha1[1].re / beta1[1].im);
          } else {
            D[5].re = alpha1[1].im / beta1[1].im;
            D[5].im = -(alpha1[1].re / beta1[1].im);
          }
        } else {
          t = std::abs(beta1[1].re);
          scale = std::abs(beta1[1].im);
          if (t > scale) {
            scale = beta1[1].im / beta1[1].re;
            absxk = beta1[1].re + scale * beta1[1].im;
            D[5].re = (alpha1[1].re + scale * alpha1[1].im) / absxk;
            D[5].im = (alpha1[1].im - scale * alpha1[1].re) / absxk;
          } else if (scale == t) {
            if (beta1[1].re > 0.0) {
              scale = 0.5;
            } else {
              scale = -0.5;
            }

            if (beta1[1].im > 0.0) {
              absxk = 0.5;
            } else {
              absxk = -0.5;
            }

            D[5].re = (alpha1[1].re * scale + alpha1[1].im * absxk) / t;
            D[5].im = (alpha1[1].im * scale - alpha1[1].re * absxk) / t;
          } else {
            scale = beta1[1].re / beta1[1].im;
            absxk = beta1[1].im + scale * beta1[1].re;
            D[5].re = (scale * alpha1[1].re + alpha1[1].im) / absxk;
            D[5].im = (scale * alpha1[1].im - alpha1[1].re) / absxk;
          }
        }

        if (beta1[2].im == 0.0) {
          if (alpha1[2].im == 0.0) {
            D[10].re = alpha1[2].re / beta1[2].re;
            D[10].im = 0.0;
          } else if (alpha1[2].re == 0.0) {
            D[10].re = 0.0;
            D[10].im = alpha1[2].im / beta1[2].re;
          } else {
            D[10].re = alpha1[2].re / beta1[2].re;
            D[10].im = alpha1[2].im / beta1[2].re;
          }
        } else if (beta1[2].re == 0.0) {
          if (alpha1[2].re == 0.0) {
            D[10].re = alpha1[2].im / beta1[2].im;
            D[10].im = 0.0;
          } else if (alpha1[2].im == 0.0) {
            D[10].re = 0.0;
            D[10].im = -(alpha1[2].re / beta1[2].im);
          } else {
            D[10].re = alpha1[2].im / beta1[2].im;
            D[10].im = -(alpha1[2].re / beta1[2].im);
          }
        } else {
          t = std::abs(beta1[2].re);
          scale = std::abs(beta1[2].im);
          if (t > scale) {
            scale = beta1[2].im / beta1[2].re;
            absxk = beta1[2].re + scale * beta1[2].im;
            D[10].re = (alpha1[2].re + scale * alpha1[2].im) / absxk;
            D[10].im = (alpha1[2].im - scale * alpha1[2].re) / absxk;
          } else if (scale == t) {
            if (beta1[2].re > 0.0) {
              scale = 0.5;
            } else {
              scale = -0.5;
            }

            if (beta1[2].im > 0.0) {
              absxk = 0.5;
            } else {
              absxk = -0.5;
            }

            D[10].re = (alpha1[2].re * scale + alpha1[2].im * absxk) / t;
            D[10].im = (alpha1[2].im * scale - alpha1[2].re * absxk) / t;
          } else {
            scale = beta1[2].re / beta1[2].im;
            absxk = beta1[2].im + scale * beta1[2].re;
            D[10].re = (scale * alpha1[2].re + alpha1[2].im) / absxk;
            D[10].im = (scale * alpha1[2].im - alpha1[2].re) / absxk;
          }
        }

        if (beta1[3].im == 0.0) {
          if (alpha1[3].im == 0.0) {
            D[15].re = alpha1[3].re / beta1[3].re;
            D[15].im = 0.0;
          } else if (alpha1[3].re == 0.0) {
            D[15].re = 0.0;
            D[15].im = alpha1[3].im / beta1[3].re;
          } else {
            D[15].re = alpha1[3].re / beta1[3].re;
            D[15].im = alpha1[3].im / beta1[3].re;
          }
        } else if (beta1[3].re == 0.0) {
          if (alpha1[3].re == 0.0) {
            D[15].re = alpha1[3].im / beta1[3].im;
            D[15].im = 0.0;
          } else if (alpha1[3].im == 0.0) {
            D[15].re = 0.0;
            D[15].im = -(alpha1[3].re / beta1[3].im);
          } else {
            D[15].re = alpha1[3].im / beta1[3].im;
            D[15].im = -(alpha1[3].re / beta1[3].im);
          }
        } else {
          t = std::abs(beta1[3].re);
          scale = std::abs(beta1[3].im);
          if (t > scale) {
            scale = beta1[3].im / beta1[3].re;
            absxk = beta1[3].re + scale * beta1[3].im;
            D[15].re = (alpha1[3].re + scale * alpha1[3].im) / absxk;
            D[15].im = (alpha1[3].im - scale * alpha1[3].re) / absxk;
          } else if (scale == t) {
            if (beta1[3].re > 0.0) {
              scale = 0.5;
            } else {
              scale = -0.5;
            }

            if (beta1[3].im > 0.0) {
              absxk = 0.5;
            } else {
              absxk = -0.5;
            }

            D[15].re = (alpha1[3].re * scale + alpha1[3].im * absxk) / t;
            D[15].im = (alpha1[3].im * scale - alpha1[3].re * absxk) / t;
          } else {
            scale = beta1[3].re / beta1[3].im;
            absxk = beta1[3].im + scale * beta1[3].re;
            D[15].re = (scale * alpha1[3].re + alpha1[3].im) / absxk;
            D[15].im = (scale * alpha1[3].im - alpha1[3].re) / absxk;
          }
        }
      }
    }
  }
}

//
// File trailer for eig.cpp
//
// [EOF]
//
